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Abstract 

We propose a modified power method for computing the subdominant eigenvalue A2 of a matrix 
or continuous operator. While useful both deterministically and stochastically, we focus on defining 
simple Monte Carlo methods for its application. The methods presented use random walkers of 
mixed signs to represent the subdominant eigenfuction. Accordingly, the methods must cancel these 
signs properly in order to sample this eigenfunction faithfully. We present a simple procedure to 
solve this sign problem and then test our Monte Carlo methods by computing the A2 of various 
Markov chain transition matrices. As the | A2 1 of this matrix controls the rate at which Monte Carlo 
sampling relaxes to a stationary condition, its computation also enabled us to compare efficiencies 
of several Monte Carlo algorithms as applied to two quite different types of problems. We first 
computed A2 for several one and two dimensional Ising models, which have a discrete phase space, 
and compared the relative efficiencies of the Metropolis and heat-bath algorithms as a function 
of temperature and applied magnetic field. Next, we computed A2 for a model of an interacting 
gas trapped by a harmonic potential, which has a mutidimensional continuous phase space, and 
studied the efficiency of the Metropolis algorithm as a function of temperature and the maximum 
allowable step size A. Based on the A2 criterion, we found for the Ising models that small lattices 
appear to give an adequate picture of comparative efficiency and that the heat-bath algorithm is 
more efficient than the Metropolis algorithm only at low temperatures where both algorithms are 
inefficient. For the harmonic trap problem, we found that the traditional rule-of-thumb of adjusting 
A so the Metropolis acceptance rate is around 50% range is often sub-optimal. In general, as a 
function of temperature or A, A2 for this model displayed trends defining optimal efficiency that 
the acceptance ratio does not. The cases studied also suggested that Monte Carlo simulations for 
a continuum model are likely more efficient than those for a discretized version of the model. 
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I. INTRODUCTION 



When designing a Monte Carlo simulation, the computational scientist often must decide 
which of several algorithmic options is the most efficient or how to optimize a particular 
algorithm. Besides experience, few rules of guidance exist. 

For detailed balance algorithms, the class we assume, likely the most rigorous rule is 
based on the work of Peskun pp. He showed analytically that if P l and P 2 are two transition 
matrices that satisfy detailed balance and asymptote to the same limiting distribution, then 
for jumps from state j to state i, algorithm 1 is more efficient than algorithm 2 if 

4>4, for alii ^j. (1) 

While he used this relation to establish the greater efficiency of the Metropolis-Hastings 
algorithm over several other generalized Metropolis algorithms, establishing this relation on 
a case by case basis is generally difficult for large phase spaces. 

The theory of Markov chains says the magnitude of the transition matrix's subdominant 
eigenvalue A2 controls the rate at which the sampling relaxes to a stationary condition. 
Doll et al. [2j call this eigenvalue the asymptotic convergence parameter j3]. Because its 
magnitude must be less than one, its closeness to zero indicates a high degree of efficiency, 
and closeness to 1, poor efficiency. Statisticians in particular [4] have derived a number of 
upper and lower bonds for this eigenvalue. Again, making use of this rigorous information 
is sometimes difficult. 

Oral tradition says that efficient sampling occurs when the acceptance ratio is around 
50%. The acceptance ratio is the number of jumps to a different state divided by the total 
number of jumps attempted. In some sense tradition is consistent with Peskun. Because 
the transition matrix satisfies J2i P%j — 1 ; moving transition probability off of the diagonal 
generally increases the acceptance. However, the more relevant implication of Peskun's 
result is the need for jumps among areas distantly separated in phase space. Acceptance 
ratios are particularly misleading in cases where phase space separates into several relatively 
localized regions of large Boltzmann weight separated by an energy activation barrier larger 
than thermal fluctuations. In this circumstance, sampling is quasi-non-ergodically confined 
to one region, and within this region the acceptance ratio may be, or may have been adjusted 
to be, the canonical 50%. 
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A direct approach to assessing the second eigenvalue is computing it numerically. The 
size of the transition matrix (or in many cases the transition operator) restricts a direct de- 
terministic approach to small systems because of the amount of computer memory required. 
Recently, Doll et al. [2j used the deterministic approach on the Metropolis algorithm's 
transition matrix for single quadradtic, quartic, double well, and triple well oscillator poten- 
tial energies. They discretized the one-dimensional phase space and prohibited very large 
displacements to create a finite matrix representation of the asymmetric transition matrix. 
They then used standard eigensystem software to compute the subdominant eigenvalue as 
a function of temperature and the Metropolis box size. The Metropolis box size is the 
maximum size of the proposed move from the current position in phase space. 

Doll et al. found several interesting results. The results for the quadratic and quartic 
wells were similar: adjusting the box size so the acceptance rate is around 50% resulted 
in the second eigenvalue being reasonably removed from unity. For the double and triple 
wells, the competition between intra- and inter-well sampling made box size optimization 
more difficult. In some cases an acceptance rate around 50% corresponded to an unfavorable 
second eigenvalue. Most interestingly, apparent activation energies in the behavior of the 
second eigenvalue existed and depended on the box size. Further, as a function of box size, 
they found structure in the second eigenvalue reflecting the expected phase space structure. 
The acceptance ratio however featured little of the informative structure seen in the second 
eigenvalue. In short, the sampling dynamics, in particular the length scale dependence of 
the apparent activation energies, contained information about the underlying structure of 
the potential energy surface about the width of the barrier region, for example. 

Extending this type of deterministic analysis to higher dimensional phase spaces, however, 
is rapidly checked by inadequate computer memory. Accordingly, we are proposing a shift 
to a new Monte Carlo method to make this extension possible. We will illustrate the 
usefulness of this new method by comparing the second eigenvalue for several Monte Carlo 
algorithms applied to the one and two-dimensional Ising models and for the Metropolis 
algorithm applied to a gas of N interacting particles in a harmonic trap. 

Computing the subdominant eigenpair (A2, l^)) is a significant shift beyond conventional 
Monte Carlo eigenanalysis. The Monte Carlo methods based on the power method have long 
computed just the dominant eigenpair (Ai, of very large matrices. For the Monte Carlo 
transition matrices, the dominant eigenvalue is known and must be one. Further, the left 
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and right hand dominant eigenvectors are also known. If P is column stochastic, that is, 
YjiPij — 1) the right eigenvector is the limiting distribution. The left eigenvector has all 
its components being the same positive number. Our method based on a modified power 
method, uses this information about the dominant eigenpair to focus on obtaining the first 
subdominant pair. It starts with the recent very large matrix multiple eigenvalue Monte 
Carlo work of Booth and Gubernatis [SHE] and adds to it an algorithm just proposed by 
Yamamoto [9]. 

Our approach is significantly different from the one proposed by Blote and Nightingale 
[TUHT2] who adapted the variational quantum Monte Carlo method [12] for computing a 
few excited states to the calculation of subdominant eigenvalues. In their modificaton of 
this variational approach, they utilized knowledge of the dominant eigenpair, and with a 
multi-parameter trial wavefunction, accomplished high precision estimates of the dynamical 
exponent for a Metropolis simulation of the two-dimensional Ising model at the bulk critical 
point. Important in this method is the quality of the trial wavefunction. Recently, Casey et 
al. [H] applied this method to a multi-variate Gaussian. 

In Section 2, we will introduce our method and in Section 3 define our models and discuss 
the details of our simulations. Section 4 contains our results for the Ising and harmonic trap 
models. Because Ising model simulations have no Metropolis box size, we instead varied 
the algorithm, using the standard single spin-reversal Metropolis method, plus single- and 
multi-site heat bath algorithms. In the last section, Section 5, we summarize our results and 
discuss future work. 

II. BACKGROUND 

Most commonly used Monte Carlo eigenvalue methods are based on the power method. 
The power method projects some starting state to the eigenpair of some matrix or operator 
A associated with the eigenvalue of largest absolute value. For a Markov chain transition 
matrix P, this eigenvalue must always be real, positive, and unity. With an initial state 
the power method iterates 



1. \<P) = A\ij) 




(2) 
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until some convergence criterion is met. Upon convergence, the eigenstate of the dominant 
eigenpair is \<ft) and the eigenvalue is \\<f)\\- Any norm may be used. In this paper, we will 
use the infinity norm, ||0|| = max, \<fi i \ 1 where <pi are the components of \<fi) in some basis. 

To find two eigenpairs, we need two starting states \<p') and | <f>"). If we were to apply 
the power method to them independently, each would independently converge to the same 
dominant eigenpair. To couple the iteration, we modify the first step of the power method 
to be 

w) = m+vW), (3) 

with the intent of converging to the dominate state j^i) and to the subdominant 
state \ip2)- 

To fix values for the 7/'s, we start with the matrix-vector form of the defining equation 
for an eigenpair 

AVi = EAiV»i (4) 
j 

and note that for any ^ ^ an exact estimator for the eigenvalue is 

* = ^ < 5 > 

and related exact estimators exist for sums of these components grouped in any number of 
overlapping or non-overlapping ways 

This observation puts a constraint on the allowed values of rj for which the sum ipi = (p^ + rj^ 
is an eigenvector. If we substitute, ipi into ^ and cross-muliply the equalities involving Ri 
and R 2 , a quadratic equation for r\ results. One root of this equation, rji, makes the eigenvalue 
estimate Ai associated with <// + r](f>'- larger than the other. As shown in [7J E] , these choices 
converge l^') and \tp") to and IV^)- It is straight-forward to generalize this method to 
compute more that two eigenpairs. 

Following Yamamoto [9] , we modified this method to converge to the subdominant eigen- 
pair using knowledge of the first. We define 

P\= E 0i A = E 0i 

P[ = E E ■\,jO J l" 2 EE XjO; r 



After a sufficient number of iterations 



W) ~ + a 2 |0 2 ) ~ Oll^l) +021-02) (8) 

|V>") ~ OiAi|0i) + a 2 A 2 |0 2 ) ~ aiAil^i) + a 2 X2\^2)- (9) 



Accordingly, we can write 



Pi ~ dial + a 2 /3! P[ ~ g^o^Ai + a 2/ 9iA 2 
P 2 ~ a x a 2 + a 2 /3 2 P2 ~ a i a 2Ai + a 2 /3 2 A 2 



(10) 



where 



«1,2 = E 01,i /5l,2 = E 02,i- (11) 



We can now solve (10) for 



and 



2 aiP 2 - a 2 Pi 
A 2 P - P' 



"i (A 2 - Ai)' 

The P's and a's are easily computed sums. In computing ot\ we use our exact knowledge 
of l^i), and in computing P', we use the exact value of Ai which is unity. The new power 
method iterates 

1. \l>)=A\<j>) 

2. Calculate ai and A 2 

3. |v)^-|V)-^i|^i) 

4. \4>) = \1>)/\M (12) 

Yamamoto shows that his procedure converges to A 2 and |-0 2 ) provided 

(Ai-A 2 )/A 1 <r 7 <(A 1 + A 2 )/A 1 . (13) 

In what follows, we chose rj to be close to the lower bound. 

Both deterministic and Monte Carlo use of this modified power method are possible. In 
the following section, we will give more details of our Monte Carlo implementation. In this 
implementation, we use the left dominant eigenvector of A, that is, the one with uniform 
positive components. This choice trivializes the computation of several sums. 
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III. MODELS AND METHODS 



A. Ising Models 

We considered both the L site one-dimensional model 

L 

E = -J2(JsiS i+1 + H Si ) (14) 
i=i 

and the L x L site two-dimensional model 

L L 

E = - J2 J2(si,jS i+1J + SijSi d+1 + Hsij) (15) 

i=l j=l 

in an external field H. The Ising variables Sj = ±1. We assumed periodic boundary 
conditions; that is, Si + i = Sj in one dimension, and Sj+L,j = Sjj+i, = Sjj = Sj+Lj+L in two 
dimensions. 

For these models, we computed the second eigenvalue for the transition matrices of multi- 
ple Monte Carlo algorithms: the single-site, spin-reversal Metropolis algorithm, a single-site 
heat bath algorithm, and the two and three-site heat bath algorithms. The algorithms are 
defined by the transition probability matrix Ps>s, which in turn defines the probability of 
jumping from state \S) to \S'). With iV equal to L or L 2 , our state is 

\S) = \s 1 ,s 2 ,...,s N ). (16) 

The jump to state \S') produced by the Metropolis algorithm has a proposed flip of the Ising 
spin at one site accepted or rejected according to 

P s >s = min [l,exp(-E(S')/kT)/ exp(-E(S)/kT)] . (17) 

The single-site heat bath algorithm transitions the state to itself or to one with the spin 
reversed at one site. If E(S) is the energy of the state with the single spin reversed and 
all remaining spins fixed and Z = exp(—E(S)/kT) + exp(—E(S)/kT), then the non-zero 
elements of P are 

P ss = exp(-E(S)/kT)/Z 

Pss = exp(-E(S)/kT)/Z. (18) 

A multiple-site heat-bath algorithm is a natural extension of the above. The single-site 
heat-bath algorithm samples one of the two spin states from the conditional Boltzmann 
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distribution of a single spin with all the other spin values fixed. A multiple-site heat-bath 
algorithm samples the state of several neighboring spins with the rest fixed. For an n-site 
algorithm one of 2 n states is selected. 

For each of these algorithms, we computed the second eigenvalue deterministically using 



standard eigensystem software for small lattices, deterministically using (12) for slightly 



larger lattices, and stochastically using (12) for still larger lattices. We will now give the 
details of the Monte Carlo approach. 

The Monte Carlo method implements the modified power method by using a collection 
of M random walkers, each specified by a weight ws and a spin configuration state \S). The 
number of these states is 2^. We represented a spin configuration by an integer S in the 
range to 2 Ar_1 where each bit of this integer corresponds to a lattice site and a plus Ising 
spin maps to a set bit and a minus spin maps to an unset one. The weight represents the 



component of the subdominant eigenstate in the spin-configuration basis (16). Monte Carlo 
becomes necessary when this number is too large for a deterministic calculation. The Monte 
Carlo method is most powerful when M <C 2 N . 

The algorithm estimates the eigenvalue by using the walkers in the exact estimators (J5]) 
and ^ not by using them in a variational estimator [HI [12]. This use needs two regions. 
In general, the choice of regions is not critical [15], as long as they are populated by a 
sufficient number of walkers so that the walkers are representative of the eigenstate in that 
region. Regions therefore do not necessarily have to be exclusive and can overlap. For the 
zero-field Ising model, our regions were all states with an up-spin majority and all states 
with a down-spin majority. We note this means states with equal numbers of up and down 
spins do not contribute to the estimation. This choice accommodated the fact that as the 
lattice size increases most of the walkers at low temperature are either in the all spins-up 
or all spins-down state. In the non-zero-field simulations, one region was the one state with 
all spins aligned with the magnetic field and several adjacent nearly all-aligned states (the 
number of adjacent states that needed to be included with the all- aligned state scaled with 
system size); the other region was all other states. 

In contrast to the dominant eigenstate, which must have only non-negative components, 
the subdominant eigenstate must have some negative components. While not essential, it is 
helpful if the initial walker weights mix plus and minus signs. While use of starting states 
whose variational energy is a significantly better approximate to the answer is possible, their 
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use was unnecessary for the present study. 

For each walker in state \S), we need to sample P$'s to produce a walker in the state \S'). 
To do this we used two lists, one for the current walkers and one for the new ones produced 
by a Monte Carlo method that samples a state \S') from the cumulative probability function 
of Ps's, that is, from S|»=o Ps"S- With this procedure we produce one new walker for every 
old one. 

The third step in the algorithm, which updates the discrete components ips of after 
adding or removing a contribution from the corresponding component of the known dominant 
eigenstate, can mix oppositely signed walkers contributing to the same state \S). It is 
essential for a correct solution to have oppositely-signed walkers cancel properly. Instead of 
using the weight cancellation methods developed in [5j El [16] , we tried a simpler method: 
After the generation of the new list of walkers is complete, we scan it, identify all walkers 
in the same state, that is, identified by the same positive integer, and replace them with 
one walker whose weight is the sum of contributing walkers' weights. This list compression 
procedure worked well. 

As the iteration progresses, the repeated matrix-vector multiplication creates some very 
large-weighted walkers and some very small-weighted ones. It is inefficient to process the 
small- weighted ones. Accordingly, after we compressed our list, we eliminated the small- 
weighted ones by a weight cut-off procedure: We scan the list, and if a \ws\ fell below e, 
then we would draw a random number £. If £ was larger than \ws\/e, we would keep this 
walker but would increase its weight to tos/e. Otherwise, we would remove it from the 
list. We used e = 1CT 3 . The weight cut-off procedure reduces the number of walkers. If 
reduced too much, say by a half, we replenished the walker population size by replacing the 
largest-weighted walkers with m = integer (ws + £) walkers with weight ws/m until the list 
size was approximately restored to its original size. 

B. Harmonic Trap 

For models in the continuum, we considered a gas of N interacting classical particles in a 
harmonic trap. We intend this to be a simple model of a "cold atom" system. The potential 
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energy is 



V T rap(X) = ^Y1 Kx i + Yl V ( X i~ X i) 




(19) 



J(#o=l 



Here, X denotes a position in phase space, that is, X = (x±,X2, ■ ■ ■ ,xn). The Xi are the 
particle displacements from the trap center. We chose v(xi — Xj) = C(\xi — Xj \ — d) 2 , a rather 
artificial interaction but one staging a competition between the tendency of particles to roam 
freely and independently within the externally-generated harmonic trapping potential and 
the tendency to position themselves within a distance of d of each other to accommodate 
their mutual interactions. We always chose d—1. 

For this model, we computed the second eigenvalue for the Metropolis algorithm. Here, 
the transition probability for jumping from position X to X' is 



The Metropolis algorithm proposes a phase space position change one particle at a time, 
say changing Xi to x\. The proposal samples an x\ from T(X',X) = H^L 1 t(x' i: Xi) where 
t(x' i: Xi) = unless = |xj + A| < T, where T is a cutoff appropriate to the potential. The 
random number £ is chosen uniformly in the interval (—A, A). The parameter A is called the 
Metropolis box size. The proposed change is accepted if exp(— V(X') / kT) / exp(— V(X) / kT) 
is greater than another random number £ chosen uniformly in [0, 1]. The parameters V and A 
may be varied widely. Their ratio is what principally affects the second eigenvalues obtained. 

Our Monte Carlo strategy was analogous to the one we used for the Ising model, but 
some details were adjusted to move from discrete to continuous states. In short, we mapped 
the continuum problem onto a discrete one. We represented each walker by a weight wx and 
phase space position X, divided phase space into cells, and formed cell lists grouping walkers 
into the cells. We then defined components for the various \ip) and \4>) states by mapping a 
combination of cell numbers and particle numbers onto a positive integer. We thus defined 
discrete states, but ones typically containing many walkers. Finally, the weights of the 
walkers in each state were replaced by their average weight. Except where otherwise noted, 
the cell width in each phase space dimension was 0.05. We observed that we could reduce 
the needed number of walkers by varying the widths. The number of walkers expected to 
occupy cells at the extreme ends of the trap, for example, is much smaller than the number 
expected to occupy cells at the bottom of the trap. Therefore, enlarging the cell widths above 



P(X', X) = T(X', X) min[l, exp(-V(X')/kT)/ exp(-V(X)/kT)]. 



(20) 
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Reduced Temperature, T*=kT/J 

FIG. 1: Comparison of the deterministic and Monte Carlo calculations of the second eigenvalues of 
the Metropolis and single-site heat bath transition matrices for a zero filed one-dimensional Ising 
model on a 10 site lattice. 

0.05 units at the far ends of the trap enabled a reduction in the number of walkers. The 
optimal size of the cells at the far ends of the trap varied with the temperature, but could 
be optimized by increasing the size of these cells to the maximum size before convergence 
could no longer be achieved. We found these sizes ranged from 0.05 to 0.5 times the core 
cell width. The infinity norm of the eigenvector needed in step 5 of the algorithm was now 
taken to be the largest absolute value of the cell weights. 

We also need to define regions to calculate the parameters needed in the updating step 
3. Our definitions were simple: The first region included all walkers with all Xi < 0, while 
the second region included all walkers with all Xi > 0. Other choices of regions are possible, 
but were not thoroughly explored. Our weight control and population resizing procedures 
were the same as those used in the Ising simulations. Assigning each walker in the same cell 
their average weight was our only weight cancellation procedure. 

IV. RESULTS 
A. Ising Models 

Figure [T] is representative of the excellent fidelity of our Monte Carlo eigenvalue pre- 
dictions. Here, as a function of the reduced temperature, we compare the Metropolis's 
and single-site heat bath's algorithms subdominant eigenvalue for a zero-field, 10-site lat- 
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FIG. 2: Second eigenvalues for the Metropolis (left) and single-site heat bath (right) transition 
matrices for one-dimensional Ising lattices in zero magnetic field computed deterministically and 
by the proposed Monte Carlo method. 



tice, computed by both deterministic (using standard eigenvalue software) and Monte Carlo 
approaches. The results for other lattice sizes and for two-site and three-site heat-bath al- 
gorithms display similar accuracy. All the deterministic calculations we present were done 



by using standard software. Computing the eigenvalue by (12) gave the same result 



In Fig. [2^i, we show the second eigenvalue for the Metropolis algorithm computed for short 




12 3 4 

Reduced Temperature, T*=kT/J 



0.6 £ 



FIG. 3: Comparisons of the deterministic calculations of second eigenvalues and of the accep- 
tance ratios for the Metropolis, single-site heat bath, and two-site heat bath algorithms for a 
one-dimensional Ising 10-site lattice in zero magnetic field. Red is for the Metropolis algorithm; 
blue, for the single-site heat bath; and green, the two-site heat bath. 
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Reduced Temperature, T*=kT/J Reduced Temperature, T*=kT/J 

FIG. 4: Deterministic calculations of the second eigenvalues of the Metropolis (left) and heat bath 
(right) algorithms for the 10-site Ising lattice as a function of reduced temperature and magnetic 
field. Red is for H/J = 0; blue, for H/J = 1; and green, for H/J = 4. 

lattices deterministically and longer lattices stochastically, and in Fig. [2]d we show the same 
for the single-site heat-bath algorithm. As in Fig. [TJ the Metropolis's algorithm eigenvalue for 
reduced temperatures greater than 1 is always lower than that of the heat-bath algorithm for 
a given lattice size. For a given length and any algorithm, the magnitudes of the eigenvalues 
always decrease monotonically with increasing temperature. For reduced temperatures less 
than one, the eigenvalues closely approach unity. The Metropolis algorithm approaches 
unity faster than the heat-bath algorithm. While not shown, we remark that the second 
eigenvalue of all algorithms increase with lattice size, pushing the approach to unity to higher 
reduced temperatures. The difference in the temperature dependences of the eigenvalues and 
acceptance ratios for the Metropolis, singe-site heat bath, and two-site heat bath algorithms 
is shown in Fig. |3j We note the high temperature crossover in the behaviors of the Metropolis 
and heat bath algorithms. 

When the external field becomes non-zero, the high temperature lattice size trends are 
similar to the high temperature trends of the zero-field case. The low temperature behavior 
however changes. In Fig. [4] we see that for both the Metropolis and heat bath algorithms 
the eigenvalue no longer approaches unity for sufficiently large fields at low temperatures. 
Further, the magnitude of the eigenvalue loses its monotonic decline as a function of tem- 
perature. 

In Fig. [5] we show results for a 2 x 2, 3 x 3, and 4x4 zero-field Ising model for the 
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FIG. 5: The behavior of the second eigenvalues and acceptance ratio of the transition matrix for 
Metropolis algorithm as a function of the reduced temperature when applied to a 2 x 2, 3 x 3, and 
4x4 Ising models in a zero magnetic field. The markers record the results for the eigenvalue, and 
the solid line, for the acceptance ratio. 

Metropolis algorithm. The markers are the eigenvalue results, and the solid line, those for 
the acceptance ratios. The behavior of the second eigenvalue here is quite consistent with 
the behavior observed in the one-dimesional case. As the reduced temperature is lowered, 
the second eigenvalue approaches unity. As the system size increases, the approach occurs at 
higher and higher reduced temperatures. In the thermodynamic limit, the two-dimensional 
Ising model has a well known critical temperature at T* = 2.2692. For the 4x4 lattice the 
second eigenvalue is effectively unity well above this temperature. The loss of efficiency in 
sampling for Ising models by the Metropolis (and heat-bath) algorithms as temperature is 
lowered is more a consequence of the increasing inefficiency of the algorithms than that of 
approaching a critical point. Results for the non-zero field case are qualitatively similar to 
those for the non-zero field one- dimensional Ising model. 

In the Ising simulations, the accuracy of our results is principally controlled by the number 
of walkers used. Typical numbers ranged from 10,000 to 50,000, as the number of sites were 
varied from 9 to 25. The number of walkers was held constant over the range of temperatures 
simulated for a given lattice size and field. 
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FIG. 6: The scaling of the second eigenvalue for the Metropolis transition matrix as a function of 
A 2 /T for N = 1 (left) and 3 (right) non-interacting particles in the harmonic trap. K = 1. 

B. Harmonic Trap 

We started by considering one particle in the harmonic trap to benchmark our basic Monte 
Carlo procedures against the deterministic single harmonic oscillator results of Doll et al. 
[2], an equivalent problem approached with the use of standard eigenvalue software. In terms 



of (19), we took N = 1 and (7 = 0. We discretized the one- dimensional space, following the 
prescription of Doll et al, selected a Metropolis box size and potential-dependent cutoffs, 
and computed the transition matrix Px',x- Both the distance S proposed for the walker to 
move and the Metropolis box size A, were discrete multiples of an underlying cell width, 
which was 0.05. In units of the cell width our cut-off distance T = 400, which was more 
than sufficient for most temperatures. We again diagonalized the transition matrix using 
standard eigensystem software and found excellent agreement between our Monte Carlo 
determination, performed with the same Monte Carlo techniques used for the Ising model, 
and our deterministic results. 

To test our continuum Monte Carlo method, we exploited the observation of Doll et 
al., proven in our Appendix, that for a power-law potential x n , all the eigenvalues of the 
Metropolis Px',x are a function solely of A n /T. As illustrated in Fig. |6j we achieve excellent 
scaling over a wide range of reduced temperatures with N = 1 ( the Doll case) and N = 3 
with C = and K — 1. We note a characteristic feature of the second eigenvalue of the 
single oscillator is a minimum value. In going from the discretized case to the continuum, 
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this value decreases. As we can see from Fig. [6] adding more non- interacting particles to the 
trap increases the continuum value. 

In Fig. [7] we present a coarse summary of the behavior of the second eigenvalue and accep- 
tance ratio as a function of box size for a two-particle trap at a low and high temperature 
for several values of the trap curvature K and various values of the coupling C between 
particles. In all cases, the acceptance ratio uniformly decreases as a function of box size. 
Its sensitivity to the coupling C depends on K. It increases as C increases, less so for the 
larger values of K . At T*, the sensitivity of the second eigenvalue on C anti-correlates with 
the sensitivity of the acceptance ratio, showing more sensitivity at the smaller values of K 
and decreasing with increasing C . At T* = 2, many of the eigenvalues trends are the same 
as those for T* . The key difference is the absence of a minimum value of the eigenvalue for 
K = 0.5 but its present at the other two K values and its location shifting toward smaller 
box sizes as K increases. A minimum occurs for K = 0.5 but is located at a box size larger 
than we simulated. 

Our final results are shown in Fig. [8} Here, we have N = 3 and a single value of the 
curvature, K = 2. The results show the same trends as the N = 2 and K = 2 case in Fig. [7} 
The minimum in the eigenvalue however shifted towards smaller box size. 

Figures [6]j8] suggest that at high temperature, the box size, which we recall is measured 
in units of the cell width, should be as large as possible. This seems reasonable [17]. We 
note the acceptance ratio will then be below the 40%. At low temperature, there is typically 
an optimal box size. For the cases presented, the acceptance ratio is around 40%. 

Besides the number of walkers, the accuracy of the results also depends on the cell width. 
The value of 0.05 was found by experimentation to be convenient both with respect to 
accuracy and efficiency Typical numbers of walkers ranged from 50,000 to 200,000 as the 
number of particles was increased from one to three. Fewer walkers could be used for higher 
temperatures than for lower ones for the same system. 

V. CONCLUDING REMARKS 

We proposed and benchmarked a numerical method for computing the subdominant 
eigenvalue A2 of a matrix or continuous operator. Based on the work of [8] and [9], this 
method can be implemented deterministically and stochastically, computes just this one 
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FIG. 7: The second eigenvalues and acceptance ratios for the Metropolis algorithm as a function 
of Metropolis box size A for two particles in the harmoninc trap for various values of the coupling 
constant C. Down the left column T* = 2; down the right, T* = 10. Across the top row K = 0.5; 
the middle row, K = 1; and the bottom row, K = 2. 
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FIG. 8: The second eigenvalue for the Metropolis transition matrix for 3 particles in the harmonic 
trap for various values of their coupling constant C as a function of Metropolis box size. K = 2 
On the left, T* = 2, while on the right T* = 10. 

eigenvalue, and requires knowing the dominant eigenpair. For Markov chain transition 
matrices, this pair is known analytically so we used the method to compute the A2 of 
various transition matrices. For such matrices < | A2 1 < 1, with small | A2 1 implying 
large Monte Carlo efficiency Specifically, we computed the A2 of the transition matrices 
for several one and two dimensional Ising models, which have a discrete phase space, and 
compared the relative efficiencies of the Metropolis and heat-bath algorithms as a function 
of temperature and applied magnetic field. Based on the A2 criterion, we found that small 
lattices appear to give an adequate picture of comparative efficiency and that the heat-bath 
algorithm is more efficient than the Metropolis algorithm only at low temperatures where 
both algorithms are inefficient. We also computed the A2 of the transition matrix of a 
model of an interacting gas trapped by a harmonic potential, which has a mutidimensional 
continuous phase space, and studied the efficiency of the Metropolis algorithm as a function 
of temperature and the maximum allowable step size A. We found that the traditional 
rule-of-thumb of adjusting A so the Metropolis acceptance rate is around 50% range is often 
sub-optimal. In general, as a function of temperature or A, A 2 for this model displays 
trends defining optimal efficiency that the acceptance ratio does not. The cases studied also 
suggest that Monte Carlo simulations for a continuum model in the continuum are likely 
more efficienct than those for a discretized version of the model. 
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Many of our results and conclusions, of course, could be limited to these specific cases; 
however, we suggest establishing their degree generality would be advantageous. For exam- 
ple, we quantified situations where focusing on acceptance ratios is ill advised. From our 
experience, this result is not very surprising. Clearly, computing A2 is more justified and 
gave more specific information about which of several algorithmic approaches or adjustments 
is likely preferable. More interestingly, if we can always capture efficiency trends for small 
systems sizes, this is a significant simplification. 

Our Monte Carlo techniques were simple. In treating more complicated problems, we 
might need some of the methods used in [HI US] plus others. Our techniques and approach 
were also quite different than those by [T01 - H2] . For example, our starting functions were 
considerably simpler to construct, and we did not need to use these functions as an im- 
portance function guiding our random walkers. Additionally, we used exact estimators for 
the eigenvalue estimates instead of variational ones. On other hand, we must be concerned 
with the cancellation of positively and negatively signed walkers. Our procedures for doing 
so however where quite simple and effective, and our success means we solved a type of 
sign problem. We do not expect our eigenvalues to be as precise as is possible with this 
other method. In most cases, our eigenvalues estimate errors were smaller than marker sizes. 
For present high precision in these estimates is unnecessary. More advantageous starting 
functions could easily be used if needed. 

The most substantive issue for future work is more fully understanding the utility of 
A2 as a metric for comparative algorithmic efficiency. Focusing on this quantity has the 
advantage of it being well defined and backed by some rigorous results. From a practical 
point of view, the most efficient algorithm is the one that produces a measurement of a 
required accuracy with the less amount of computer time. Using A2 does not address the 
time is takes to perform a Monte Carlo step; it only says something about the number of 
steps needed. Experience has shown that the fast or slow relaxation of an algorithm is 
generally accompanied by fast of slow generation of statistically independent measurements. 
Likely a good A2 is a necessary but not sufficient rule-of-thumb. We note that Peskun's 
result ([T]) insures the variance of any measurable of algorithm 1 will be less than of equal to 
that computed from algorithm 2. The connection of A2 to his result about measurables is 
indirect. 
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Appendix: Scaling 



The basic equation for a Markov chain is 

p(x) = / dyP(x,y)p(y) 



Here, p is the stationary probability density of the chain and P is the transition probability 
density defining the chain. They satisfy 

J dxp(x) =1 
J dx P (x, y) = 1 

The Metropolis algorithm makes a specific choice for the transition probability density 

P (x, y) —T (x, y) min {1, exp [-0V (y)]/exp [-0V (x)]} 

where / dxT(x,y) = 1. For the proposal density T(x,y) the Metropolis choice has it non- 
zero only on the interval y — ^, y + ^ over which its magnitude is 1/A. We will show that 
for potentials of the form x n the product of this choice of the proposal probability density 
and the standard Metropolis acceptance density 

A(x, y) = min {1, exp [-0V (y)]/exp [-0V (x)}} 

lead to a scaling of P(x,y) that implies all it eigenvalues scale as a function of A n /T 

We will measure all lengths in units of A and start by writing the acceptance function as 
A(x,y) = A(x,y;T) to make the T dependence explicit. Next, for potentials of the type x n 



A(x,y;T) = min jl, exp (-|^^/exp 



x 



min < 1, exp 



A I- V -- T - 
\A' A' A" 



A" J 



exp 



(a«) 



Thus, from the point of view of the acceptance, going from y to x at temperature T is the 
same as going from x/A to y/A at T/A n . 

We now write T (x, y) = T (x, y; A) which distributes x uniformly over an interval cen- 
tered at y and of width A, that is, 

f A A 

xe^y- — <y<y+ — 
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Over this interval its amplitude is 1/A. Next, we consider 

x_ fy _i< y_ < y_ U 
A LA 2 - A - A 2i 

A function that distributes x/A uniformly over this unit interval and has a unit amplitude 
will make an acceptable proposal probability. Such a function is A T (j^, ^; . Thus for 
the scaled system we have 

p ( — —■ 1 = A T f — V--i\a(- 

VA'A' 'AW ' U'A' j U'A'A-j 

= T(x,y;A)A(z,y;T) 
= P(x,y;A,T) 

which establishes the scaling of all the eigenvalues of P (x, y; A, T) 
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